Diffusion Monte Carlo study of circular quantum dots 
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We present ground and excited state energies obtained from Diffusion Monte Carlo (DMC) 
calculations, using accurate multiconfiguration wave functions, for N electrons (JV < 13) confined 
to a circular quantum dot. We compare the density and correlation energies to the predictions 
of local spin density approximation theory (LSDA), and Hartree-Fock theory (HF), and analyze 
the electron-electron pair-correlation functions. The DMC estimated change in electrochemical 
potential as a function of the number of electrons in the dot is compared to that from LSDA and 
HF calculations. Hund's first rule is found to be satisfied for all dots except N = 4 for which there 
is a near degeneracy. 
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I. INTRODUCTION 



Modern microfabsjeation technology is capable of mak- 



as* 



ing quantum dotaH'H that are sufficiently small that 
they contain only a small number of mobile electrons. 
There has been much interest in studying the atomic^, 
like properties q£ these dots with tunnel conductance^ 
and capacitances experiments. The ground states of 
clean circular dots exhibit, shell structure and believed 
to obey Hund's first rulefltil. The shell structure is par- 
ticularly evident in measurements of the change in elec- 
trochemical potential due to the addition of one extra 
electron A^v = fi(N + 1) — fi(N) where N is the number 
of electrons in the dot, and fJ,(N) = E(N) - E(N - 1) is 
the electrochemical potential of the system. Theoretical 
predictions of Ajy and the excitation energy spectrum 
require accurate calculations ground-state and excited- 
state energies. Exact diagonalization studieaffl are accu- 
rate for a very small number of electrons, but the num- 
ber of basis functions needed to obtain a given accuracy 
and the computational cost grow very rapidly with elec- 
tron number.|-.Jn practice they have been used for up 
to 8 electronsQa, but thg accuracy is very limited for all 
except N < 3. Hartreecl, restricted Hartre&^ock (HF), 
spin-and/or-space unrestricted Hartree-FockEjES (UHF) 
and local spin-density (LSDA) and current density func- 
tional methodsEjtJ, give results that are satisfactory for 
a qualitative understanding of some systematic proper- 
ties. However, comparisons with exact results show dis- 
crepancies in the energies that are substantial on the 
scale of energy differences. An advantage of the approx- 
imate approaches is that no serious size and geometry 
constraints are imposed. 

In this paper we employ the Quantum Monte Carlo 
(QMC) method (both variational Monte Carlo (VMC) 



and diffusion Monte Carlo (DMC)) because they yield 
very accurate energies at a computational cost that grows 
relatively modestly with the number of electrons. The 
statistical error of these calculations can be made small, 
even for dots with several tens of electrons, within a rea- 
sonable amount of computer time on a modern worksta- 
tion. In addition to the statistical error there is a sys- 
tematic error due to using the fixed-node approximation. 
This error can be reduced by optimizing the trial wave- 
functions. For the trial wavefunctions used in the present 
work, the fixed-node errors are small compared to the er- 
rors of other approximate methods. This is demonstrated 
by performing internal checks within the method and by 
comparing to the fe^-energies, available from exact di- 
agonalization studiesoH for small dots, that are accurate 
enough to make a meaningful comparison. Hence our 
results can be regarded as a benchmark to assess the ac- 
curacy of other approximate methods. In particular we 
find that, in contrast to the situation with atoms, the 
energies obtained from the LSDA method are consider- 
ably more accurate than those from the HF method. The 
same is true for the spin-densities in those cases where 
the LSDA wavefunctions are eigenstates of the total spin 
operator S 2 . 

Earlier QMC calculations on quantum dots include 
VMC calculations for circula^-dotslHI and DMC calcu- 
lations for 3-dimensional dotsEa. Fixed-phase DMC has 
been applied to dotst2l with N < 4. Path-integral Monte 
Carlo calculations have been performedo for dots with 
N < 8 but the results of these calculations bear no re- 
semblance ^either our results or those from exact diag- 
onalizationQtl. 
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II. COMPUTATIONAL METHOD 

A. Hamiltonian 

The usual model! for a disk-shaped vertical quantum 
dot is a 2-dimensional system of N electrons moving in 
the z = plane, confined by a parabolic lateral confining 
potential V con (r). The Hamiltonian is 

"'B-^ + Wr,))^ . (1) 

i— 1 i<j ' J ' 

In Eq. (0), to* is the electron effective mass, and e is 
the dielectric constant of the semiconductor. In the fol- 
lowing (if not explicitly specified otherwise) we will use 
effective atomic units, defined by h = e 2 /e = m* = 1. In 
this system of units, the length unit is the Bohr radius ao 
times em e /m*, and the energy unit is the Hartree times 
in* j (m e e 2 ). For the GaAs dots we consider here, e = 12.4 
and to* = 0.067TO e , and the effective Bohr radius a*, and 
effective Hartree H* are ~ 97.93A and ~ 11.86 meV re- 
spectively. In this first application of the method, we will 
consider circular dots with N < 13, and a parabolic po- 
tential, V con (r) = m*u 2 r 2 /2 (huj=0.28 H* = 3.32 meV), 
which should approximate the experimental situation in 
Ref. H Extensions of the calculation to N > 13, mag- 
netic field B ^ and a non-parabolic confining potential 
are in progress. 

Comparison of energies and other quantities with those 
in the literature are complicated by the fact that various 
authors use different values for the parameters, m*,e,w 
in the Hamiltonian. Note however, that two Hamilto- 
nians, Hi and H2, characterized by to*,wi and t\ and 
to*, , Ct>2 and £2 , respectively, must have the same energy 
spectrum aside from a multiplicative scale factor, i.e., 
Eu/E2i = m^/ml — lui/u>2 — £2/^1, where i labels 
the energy states of a given Hamiltonian. An interest- 
ing aspect of quantum dots is that it is possible to tune, 
A, the dimensionless ratio of the Coulomb interaction 
strength to the confining potential A = [e 2 /(elo)]/hu, 
where Iq — \f%] (m*uj), thereby allowing one to study 
both weakly interacting and strongly interacting cases. 
Our present calculations are for A = 1.89. 

B. Quantum Monte Carlo Methods 

One advantage of the QMC methods is that no restric- 
tion is placed on the form of the trial wavefunction. In 
the VMC method, Monte Carlo integration is used to cal- 
culate the many-dimensional integrals, and the parame- 
ters in the trial wavefunction can be varied to minimize 
the energy or the fluctuations of the local energy. More 
accurate results can be obtained from the fixed-node dif- 
fusion Monte Carlo method which projects, from an an- 
tisymmetric trial wave function, the lowest energy state, 
consistent with the boundary condition of preserving the 



same nodal surface. In the limit that the trial wave- 
function has the correct nodes, fixed-node DMC yields 
the exact energy with only a statistical error that can 
be made arbitrarily small by increasing the number of 
Monte Carlo steps. A detailed description of our imple- 
mentation can be found in Ref. ^l[ The fixed-node error 
is usually small compared to the errors from other meth- 
ods but it is unknown except in those cases where exact 
results are available. 



C. Trial Wavefunctions 

The errors of VMC and fixed-node DMC calculations 
depend on the quality of the trial wavefunctions. The 
trial wavefunctions we use have the form 

*(R)l,s = exp[0(R)] ' S ( R ) . (2) 

i=l 

where R = {ri ■ • • r?v} are the coordinates of the N elec- 
trons in the dot, and the «i are variational parameters. 
The configuration-state-functions, S L,S , are eigenstates 
of the total angular momentum L = L z with eigenvalue 
L and of the total spin S 2 with eigenvalue S(S + 1), and 
have the following form: 

Sf' S = £/V?jl?j , (3) 
i=i 

where the are Slater determinants of spin-up and 
spin-down electrons, using orbitals from a local density 
approximation (LDA) calculation with the same confin- 
ing potential and the same number of electrons. The 
nii are the number of determinants in the i— th configu- 
ration. In general the D^Dj are not eigenstates of S 2 . 
The coefficients (3 in the linear combination of Eq. (3) 
are fixed by diagonalizing S 2 in that determinantal ba- 
sis. For N < 13, the number of configurations, Ar con f 
and Slater determinants Ad t = Yli^T* m ii appearing in 
Eqs. 2 and 3 are shown in Table | and were determined 
by limiting the basis space to spin-up and spin-down or- 
bitals with |n, I) = |0, 0) for N < 2 dots, \n, I) = |0, 0) and 
|0,±l),for 3 < N < 6 dots, \n,l) = |0,0), |0,±1), |0,±2), 
and |1,0), for 7 < N < 12 dots, and. \n,l) = |0,0), 
|0,±1), |0,±2), 1 1, 0) , |0,±3), and |1,±1) for the N = 13 
dot. The noninterapting single-particle energy levels are 
e n> j = (2n+|/| + l)tJ 2 J. Basis states are then built by con- 
sidering all possible occupations of open-shell levels. For 
example, in the case of the N = 9 dot, the first 6 electrons 
fill the |0, 0) and |0, ±1) orbitals and are considered to be 
core electrons in a closed shell. Then, the wave function 
for the state \L = 0,S = 1/2) of the N — 9 dot has 3 
open-shell electrons, and includes two (iV con f = 2) con- 
figuration state functions which are linear combinations 
of toi = 2 and TO2 = 3 Slater determinants respectively. 
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The function exp[<jf>] in Eq. ( t2|) i s a generalized Jastrow 
factor of the form used in Ref. 
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2 \ l + b( ri )r i:j l + 6(rj)r,. 



(4) 



where 



6(r) = b% + b\ j tan- 1 !^ - R C ) 2 /2R C A}. 



(5) 



It explicitly includes one- and two-body correlations and 
effective multi-body correlations through the spatial de- 
pendence of b(r). The quantity R c represents an "ef- 
fective" radius of the dot, and has been assumed to be 
equal to 1.93^/N. The bo and b\ parameters depend only 
on the relative spin configuration of the pair ij. The 
parameters ay are fixed in order to satisfy the cusp con- 
ditions, that is, the condition of finitcness of the local 



energy for 



0. For a two dimensional sys- 



tem, ay = 1 if the electron pair ij has antiparallel spin, 
and ay = 1/3 otherwise. The dependence of ay on the 
relative spin orientation of the electron pair introduces 
spin-contamination into the wavefunction. However, the 
magnitude of the spin contamination and its effect on the 
energy has been shown to be totally negligible in the case 
of well optimized atomic wavefunctionsB-3 and we expect 
that to be true here as well. 

The coefficients m the one-body term, the coeffi- 
cients A, bo, and b\ in the two-body term and the coef- 
ficients (Xi multiplying the configuration state functions 
are aatimized by minimizing the variance of the local en- 
ergyEJl The resulting wavefunctions had rms fluctuations 
of the local energy that range from 0.021 H* for N=2 to 
0.255 H* for N=13. 



III. RESULTS 

Using these optimized wave functions for importance 
sampling, we perform fixed-node diffusion Monte Carlo 
calculations. We attempt to establish the accuracy of 
the fixed-node energies obtained with our trial wavefunc- 
tions by comparing them to energies from exact diago- 
nalization studies and also by performing internal checks 
within our calculations. Unfortunately, although there 
exist several papers on exact diagonalizationLrcl, the re- 
sults are usually presented in plots, rather than in tables. 
The only number we know of is in Ref. 0, who give 
an energy of 26.82 meV far N = 3, that they credit to 
Hawrylak and PfannkucheO. Starting from a single Slater 
determinant of LDA orbitals (constructed from the spin- 
up and spin-down \n, I) = 1 0, 0) states and the spin-up 
1 0, 1) spin up) state, we obtain a fixed-node DMC en- 
ergy of 26.8214(36) meV, using their model parameters 
(to* = 0.067m e ,e = 12.4 and hu = 3.37 eV), which 



is indistinguishable from the exact energy to the num- 
ber of digits quoted. We attempted also to estimate the 
fixed-node error by varying the orbitals in the determi- 
nants and by varying the number of determinants. For 
N = 7, QMC calculations using local density approxima- 
tion (LDA) and local spin-density approximation (LSDA) 
orbitals were performed. The LSDA orbitals yielded bet- 
ter VMC results (the energy was lowered by 57 mH* , and 
the fluctuations of the local energy by 7 mH*) but the 
DMC energies were unchanged within statistical uncer- 
tainty. We checked the dependence of the energy on the 
number of configuration state functions for the first ex- 
cited state of the N = 9 dot. Somewhat to our surprise, 
energies of the one-configuration (three-determinant) and 
the two-configuration (five-determinant) wavefunctions 
agreed to within 1 mH*, not only within DMC but also 
within VMC. 



A. Ground-state energies 

The ground state energies are listed in Table [n] and 
compared with results of HF, and LSDA calculations us- 
ing the Tanatar-Ceperley parametrization for the corre- 
lation energyEl The HF energies are 0.12 to 0.97 H* 
higher than the DMC energies whereas the LSDA en- 
ergies are only 0.021 to 0.042 mH* (0.25 to 0.50 meV) 
higher. In contrast, in atoms and molecules the Hartree- 
Fock total energy is considerably better than the LSDA 
total energy. There are two likely reasons for this differ- 
ence. First, Hartree-Fock treats exchange exactly while 
completely ignoring correlation, whereas in LSDA both 
exchange and correlation are approximated. In atoms 
and molecules, the exchange energy, E x , is much larger 
than the correlation energy, E c , but for the dots it is 
not, e.g. E x /E c w 30 for a Neon atom but E x /E c w 4 
for a N — 10 dot. The second reason is that the dots are 
more homogeneous than atoms or molecules and so the 
local-density approximations to E x and E c work better. 
Note also that the HF errors increase monotonically with 
electron number but the LSDA errors do not show any 
obvious trend. 



B. Excited-state energies 

In Table | we list the low-lying excitation energies for 
the N = 4... 11 dots. Koskinen et al. (Ref. |l|) find 
that the lowest excitation energies from LSDA calcula- 
tions are 11.5 mH* and 2.31 mH* for N = 8, 10 dots 
respectively, whereas our DMC calculations show that 
the corresponding lowest excitation energies are 22 mH* 
and 2 mH*. They claim that these lowest excited states 
have a spin density wave even though S = 0, but in fact 
this is just an artifact due to these LSDA wavefunctions 
not beingreigenstates of S 2 , as pointed out by Hirose and 
Wingreenli-3. For both the N = 8 and the N = 10 dots, 
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the first excited state LSDA wavefunctions are in fact lin- 
ear combinations of (L, S) = (2, 1) and the (L, S) = (2, 0) 
wavefunctions. In general, the single-determinant LSDA 
wavefunctions are eigenstates of S z , but they are eigen- 
states of S 2 only when | S z | has the maximum value con- 
sistent with filling the lowest N/2 orbitals and the exclu- 
sion principle. In other cases it is necessary to have more 
than one determinant in order to have the correct spin 
symmetry. Single-particle levels in a parabolic potential 
with the same value of 2n + \l\ + 1 are degenerate. How- 
ever, the self-consistent LSDA potential is not parabolic 
and consequently of two levels with the same 2n+ \l\ + 1, 
the LSDA orbital with the larger value of | / 1 is lower than 
the other. This serves to explain the ordering of levels for 
the dots where N differs from a closed shell by one. For 
example, the (L,S) = (2,1/2) state lies lower than the 
(L, S) = (0, 1/2) state in the N = 7 dot but the order is 
reversed in the n = 11 dot because the (n, I) — (0, ±2), 
LSDA single-particle level lies below the (n,l) = (1,1) 
level. 



C. Change in electrochemical potential 



The DMC estimates for the change in electrochemi- 
cal potential An (in meV) as a function of N are re- 
ported in Fig. [l] together with those from LSDA and HF 
calculations. We see structures and peaks at electron 
numbers 2,4,6,9 and 12 in agreement with the experi- 
ments of Ref. ^. In the independent-particle model with 
a parabolic potential, Ajv has peaks of magnitude u> at 
N = 2,6,12..., corresponding of closed shells, and is 
elsewhere. Additional features are due to the electron- 
electron interaction. It is difficult to make a more de- 
tailed comparison between experiments and theory be- 
cause of uncertainties in the Hamiltonian. In particular, 
the external potential may not be strictly parabolic and 
our assumption that u) is independent of N may not be 
an accurate description of the experimental situation. 



D. Correlation energies 

In Fig. we plot the DMC correlation energy calcu- 
lated as the-difference between the DMC energy and the 
HF energycj as a function of the electron number N. 
The dashed line indicates the LSDA correlation energy. 
From the figure one sees that LSDA overestimates the 
correlation energy by 10%-15% almost independently of 
TV. The LSDA overestimate of the correlation energy is 
smaller than in atoms and jellium spheres, where it is as 
much as 100% (Ref. fify and 30% (Ref. ||), respectively. 



E. Hund's first rule 

From Table [| we see that Hund's first rule, according 
to which the total spin of the ground state takes the max- 
imum value consistent with electrons being in the same 
shell and the exclusion principle, is satisfied for all val- 
ues of TV studied in this work, except for TV = 4. For 
TV = 4 the \L,S) = |0,0) state is just 2.2 mH* or 0.026 
meV lower than the |0, 1) state, so a small change in 
the Hamiltonian, e.g., an increase in the spring-constant 
of the confining potential, u>, could alter the ordering of 
these two states. Our result for the TV = 4 dot is in. 
qualitative agreement with the QMC results of Boltont£l 
but they find that the the singlet state is lower than the 
triplet by a larger amount (1.5 meV) than we do, for 
Hamiltonian parameters that are close to, but not equal 
to, the the ones we use. However, our result disagrees 
with our LSDA calculations which find no violations of 
Hund'sajule for TV < 13 as well as the earlier LSDA calcu- 
lationstS which found that, for dots with even TV, Hund's 
rule is satisfied for TV < 22, but violated for TV = 24. On 
the other hand spin-andi-space unrestricted Hartree Fock 
(sS-UHF) calculations!!^ predict that Hund's rule is vio- 
lated not only for TV = 4 but also for TV = 8 and TV = 9. 
It should be noted that the sS-UHF were performed for a 
smaller value, A = 1.48, of the dimensionless ratio of the 
Coulomb interacti on st rength to the confining potential, 
defined in Section II A , than our calculations which were 
for A = 1.89. Since, according to Ref. [l^, Hund's rule 
violations are less likely for smaller values of A, it is clear 
that the difference is not due to the different value of A. 
Experimental evidence indicates that Hund's rule is sat- 
isfied for TV — 4 circular dotstffl but that a small elliptical 
deformation-is sufficient for the singlet and triplet ener- 
gies to crossO, thereby confirming our finding that the two 
states are very close in energy. Given the uncertainty in 
the experimental Hamiltonian and the near degeneracy 
of the two states, it is not surprising that we find that 
the the singlet state is lowest whereas the experimental 
finding is that the triplet is lowest for the circular dot. 
The results of two exact diaganalization studies are also 
relevant in this context: EtoQ found that a small mag- 
netic field is sufficient to switch the order of the states 
whereas Hirose and WingreenO find that a small quartic 
term in the Hamiltonian has the same effect. 



F. Spin densities 

In Fig. ||we compare the spin densities, p^ and p±, and 
the magnetization, m(r) = p-\(r) — Pi(r) for the N = 9 
ground state obtained from DMC and LSDA. In this case 
the LSDA wavefunction is an eigenstate of S 2 . The agree- 
ment of the curves is impressive and extends to the whole 
region of r, including the edge, where the density gradi- 
ents are large. The same kind of agreement with LSDA 
was also obtained in the case of variational Monte Carlo 
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densititesE3 of jellium spheres. In general, it appears that 
LSDA gives accurate spin densities in those cases that the 
Kohn-Sham wavefunction has the correct spin symmetry. 
In contrast the HF spin densities show much larger oscil- 
lations than the DMC spin densities. The same behavior 
has previously been noticed for atoms, but to a much 
lesser degreecM 

Yannouleas and Landman plot in Fig. 2 of Ref. [l^ the 
charge density of a closed-shell N = 6 dot obtained from 
a spin-and-space unrestricted Hartree Fock (sS-UHF) cal- 
culation. They find that the charge density of a dot 
with a dimensionless interaction strength of A = 1.48 
has a non circular charge density that they refer to as 
a Wigner crystallized state, although the usual defini- 
tion of Wigner crystallization refers to the occurrence of 
long-range order in the 2-body density rather than short- 
range order in the 1-body density. Since the ground state 
of the N = 6 dot is of 1 S symmetry, it is apparent that 
the density must be circularly symmetric and their re- 
sult is an artifact of their computational method. In this 
context it should be noted that for very large values of A 
one cannot immediately rule out the possibility that the 
single-particle picture breaks down completely and that 
the ground state is not of 1 S symmetry. Also, in the pres- 
ence of a strong magnetic field the single-particle levels 
will reorder and the ground state need not have S sym=. 
metry. Finally, it should be noted that other authorsO 
have considered models for dots in which the confining 
potential itself can deform and therefore not be circularly 
symmetric. In this case, of course, the ground state den- 
sity of the N — 6 dot need not be circularly symmetric 
either. 



G. Pair correlation functions 

In Fig. [| we show the spherical average of the electron- 
electron pair correlation functions <7<ti,<t 2 ( r ii r 2) in the 
N = 9 case. The different behavior for pairs with parallel 
and antiparallel spin is due to the fact that the wave- 
function vanishes when parallel-spin electrons coalesce 
but not when antiparallel-spin electrons coalesce. For 
N = 9, it follows from Hund's rule that there are twice 
as many up-spin electrons as down-spin electrons. This 
is reflected in the shape of the gjj and gfj_j curves. 



IV. CONCLUSIONS 

In conclusion, we have calculated QMC ground-state 
energies, excitation energies, correlation energies, change 
in electrochemical potential due to adding an electron, 
spin densities and pair-correlation functions for circular 
quantum dots with N < 13 electrons and compared them 
to the corresponding quantities obtained from HF and 
LSDA calculations. We find that HF energies are in er- 
ror by 0.12 to 0.97 H* but LSDA energies by only 0.021 



to 0.042 H*. However, even the LSDA energies are not 
sufficiently accurate to give reliable excitation energies 
or changes in the electrochemical potential. The LSDA 
correlation energies differ from the DMC ones by~ 10%- 
15%. Hund's first rule is found to be satisfied for all dots 
up to N < 13, for the Hamiltonian parameter values em- 
ployed, except for N = 4 which has a near degeneracy. 
The LSDA spin densities are in remarkably good agree- 
ment with DMC for those cases where the Kohn-Sham 
wave function is an eigenstate of S 2 but the HF densities 
have oscillations that are too large. Finally, the pair- 
correlation functions may be of utility in constructing 
more accurate energy density functionals than the LSDA 
for 2-dimensional systems. 
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TABLE I. Ground state energies (in H*) and low- lying ex- 
citation energies (in mH*) for N < 13 dots. Also shown are 
the quantum numbers of the states and the number of config- 
uration state functions iV con f and the number of determinants 
iVdet used in constructing them. The numbers in parentheses 
are the statistical uncertainties in the last digit. 
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TABLE II. Comparison of ground state energies (in H*) 
for the dots with 2 < N < 13 computed by Hartree-Fock, 
LSDA, VMC and DMC. Also shown are the LSDA errors 
in the energy, AElsda = -Blsda - Edmc, which are much 
smaller than the HF errors Ehf — £dmc- The numbers in 
parentheses are the statistical uncertainties in the last digit. 
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FIG. 2. Correlation energies E c for circular dots computed 
with DMC (filled squares) and LSDA (open circles). 
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FIG. 1. Change in electrochemical potential A at as a func- 
tion of the number of electrons, N, in the dot. The numbers 
in the plot are the DMC spin polarizations 2S Z =N f —N J.. 
The LSDA and HF spin polarization for the N = 4 dot is 
given in parentheses; for all other N they are the same as for 
DMC. 




-0.02 



LSDA 

HF 

DMC 



10 



FIG. 3. Spin densities, pj(r), /?|(r) and magnetization 
m(r) = Pf(r) — Pi(r) as a function of distance from the cen- 
ter for the ground state of the N = 9 dot. Solid lines, DMC; 
dotted-dashed lines, LSDA, dashed line, HF. The LSDA spin 
densities for this state agree well with the DMC spin densi- 
ties but the HF spin densities have considerably larger oscil- 
lations. 
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FIG. 4. Electron pair correlation functions, from varia- 
tional MC, of the ground state of the N = 9 dot. The first 
electron is at the center of the dot. Empty circles: g(r)n; 
filled squares: g(r)if, crosses: g(r)u- 
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